COMPUTING STABILITY OF MULTI-DIMENSIONAL 
TRAVELLING WAVES* 
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Abstract. We present a numerical method for computing the pure-point spectrum associated 
with the linear stability of multi-dimensional travelling fronts to parabolic nonlinear systems. Our 
method is based on the Evans function shooting approach. Transverse to the direction of propa- 
gation we project the spectral equations onto a finite Fourier basis. This generates a large, linear, 
one-dimensional system of equations for the longitudinal Fourier coefficients. We construct the stable 
and unstable solution subspaces associated with the longitudinal far-field zero boundary conditions, 
retaining only the information required for matching, by integrating the Riccati equations associated 
with the underlying Grassmannian manifolds. The Evans function is then the matching condition 
measuring the linear dependence of the stable and unstable subspaces and thus determines eigen- 
values. As a model application, we study the stability of two-dimensional wrinkled front solutions 
to a cubic autocatalysis model system. We compare our shooting approach with the continuous 
orthogonalization method of Humpherys and Zumbrun. We then also compare these with standard 
projection methods that directly project the spectral problem onto a finite multi-dimensional basis 
satisfying the boundary conditions. 
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1. Introduction. Our goal in this paper is to present a practical extension of 
the Evans function shooting method to computing the stability of multi- dimensional 
travelling wave solutions to parabolic nonlinear systems. The problem is thus to 
solve the linear spectral equations for small perturbations of the travelling wave. 
Hence we need to determine the values of the spectral parameter for which solutions 
match longitudinal far-field zero boundary conditions and, in more than one spatial 
dimension, the transverse boundary conditions. These are the eigenvalues. There are 
two main approaches to solving such eigenvalue problems. 

• Projection: Project the spectral equations onto a finite dimensional basis, 
which by construction satisfies the boundary conditions. Then solve the re- 
sulting matrix eigenvalue problem. 

• Shooting: Start with a specific value of the spectral parameter and the correct 
boundary conditions at one end. Shoot/integrate towards the far end. Ex- 
amine how close this solution is to satisfying the boundary conditions at the 
far end. Re-adjust the spectral parameter accordingly. Repeat the shooting 
procedure until the solution matches the far boundary conditions. 

The projection method can be applied to travelling waves of any dimension. Its 
main disadvantage is that to increase accuracy or include adaptation is costly and 



*VL is a postdoctoral fellow of the Fund of Scientific Research — Flanders (F.W.O.— Vlaanderen). 
SJAM and JN were supported by EPSRC First Grant GR/S22134/01. SJAM was also supported 
by Nuffield Foundation Newly Appointed Science Lecturers Grant SCI/180/97/129 and London 
Mathematical Society Scheme 2 Grant 2611. JN was also supported by Australian Research Council 
grant DP0559083. 

^Vakgroep Toegepaste Wiskunde en Informatica, Ghent University, Krijgslaan 281-S9, B-9000 
Gent, Belgium (Veerle.Ledoux@ugent.be) 

t Maxwell Institute of Mathematical Sciences and School of Mathematical and Computer Sciences 
Heriot-Watt University, Edinburgh EH14 4AS, UK (S. J.MalhamOma.hw.ac.uk) 

§ School of Mathematics, University of Leeds, Leeds, LS2 9JT, UK (jitse@maths.leeds.ac.uk) 

^Fakultat fiir Mathematik, Universitat Bielefeld, 33501 Bielefeld, Germany (thuemmle@math.uni- 
bielefeld.de) 



1 



2 



Ledoux, Malham, Niesen and Thiimmler 



complicated as we have to re-project onto the finer or adapted basis. Further, from 
the set of eigenvalues that are produced by the resulting large algebraic eigenvalue 
problem, care must be taken to distinguish two subsets as follows. In the limit of 
vanishing discretization scale, one subset corresponds to isolated eigenvalues of finite 
multiplicity, while the other subset converges to the absolute or essential spectrum, 
depending on the boundary conditions (see Sandstede and Scheel [26]). The latter 
subset can be detected by changing the truncated domain size and observing which 
eigenvalues drift in the complex spectral parameter plane. In contrast, for the shooting 
method, it is much easier to fine-tune accuracy and adaptivity, and zeros of the 
matching condition are isolated eigenvalues of finite multiplicity (we do not have to 
worry about spurious eigenvalues as described above). Its disadvantage though is that 
by its very description it is a one-dimensional method. 

When implementing the shooting method, a Wronskian-like determinant measur- 
ing the discrepancy in the boundary matching condition in the spectral problem is 
typically used to locate eigenvalues in the complex spectral parameter plane. This is 
known as the Evans function or miss-distance function — see Alexander, Gardner and 
Jones m or Greenberg and Marietta 13J. It is a function of the spectral parameter 
whose zeros correspond to eigenvalues, i.e. a match. In practice, we match at a point 
roughly centred at the front. Further for parabolic systems it is typically analytic in 
the right-half complex parameter plane and the argument principle allows us to glob- 
ally determine the existence of zeros and therefore unstable eigenvalues by integrating 
the Evans function along suitably chosen contours. 

The one-dimensional label of the Evans function has started to be overturned and 
in particular Deng and Nii [9], Gesztesy, Latushkin and Makarov [11], and Gesztesy, 
Latushkin and Zumbrun [12] have extended the Evans function theoretically to multi- 
dimensions. Importantly from a computational perspective, Humpherys and Zum- 
brun 119j outlined a computational method crucial to its numerical evaluation in 
the multi-dimensional scenario. Their method overcame a restriction on the order 
of the system that could be investigated numerically using shooting, in particular 
using exterior product spaces. It is based on continuous orhogonalization and is 
related to integrating along the underlying Grassmannian manifold whilst evolving 
the coordinate representation for the Grassmannian. More recently Ledoux, Mal- 
ham and Thiimmler [21] provided an alternative (though related) approach to that 
of Humpherys and Zumbrun. This approach is based on chosen coordinatizations of 
the Grassmannian manifold and integrating the resulting Riccati systems. 

The main idea of this paper is to determine the stability of multi-dimensional 
travelling waves by essentially combining the two main approaches in the natural 
way: we project the spectral problem transversely and shoot longitudinally. The 
basic details are as follows. 

1. Transverse to the direction of propagation of the travelling wave we project 
onto a finite Fourier basis. This generates a large, linear, one-dimensional 
system of equations for the longitudinal Fourier coefficients. 

2. We construct the stable and unstable solution subspaces associated with 
the longitudinal far-field zero boundary conditions by integrating with given 
Grassmannian coordinatizations (chosen differently for each subspace). 

3. The Evans function is the determinant of the matrix of vectors spanning 
both subspaces and measures their linear dependence. We evaluate it at an 
intermediate longitudinal point. 
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We compare the construction of the stable and unstable solution subspaces by in- 
tegrating the spectral problem using the continuous orthogonalization method of 
Humpherys and Zumbrun. We also compare these with standard projection methods 
that project the spectral problem onto a finite multi-dimensional basis satisfying the 
boundary conditions and solve the resulting large algebraic eigenvalue problem. 
The parabolic nonlinear systems that we will consider here are of the form 

dtU ^ B AU + cd^U + F{U), (1.1) 

on the cylindrical domain M x T. Here U is a. vector of component fields and F(U) 
represents a nonlinear coupling reaction term. The diagonal matrix B encodes the 
positive constant diffusion coefficients. We suppose we have Dirichlet boundary con- 
ditions in the infinite longitudinal direction x € M and periodic boundary conditions 
in the transverse direction y e T. We assume also, that we are in a reference frame 
travelling in the longitudinal direction with velocity c. Though we assume only two 
spatial dimensions our subsequent analysis in the rest of this paper is straightfor- 
wardly extended to higher dimensional domains which are transversely periodic. We 
also suppose we have a travelling wave solution U — Uc{x,y) of velocity c satisfying 
the boundary conditions, whose stability is the object of our scrutiny. 

As a specific application we consider a cubic autocatalysis reaction-diffusion sys- 
tem that admits wrinkled cellular travelling fronts. They appear as solutions bifur- 
cating from planar travelling wave solutions as the autocatalytic diffusion parameter 
is increased. The initiation of this first bifurcation is known as the cellular instability 
and has been studied in some detail in the reaction kinetics and combustion literature; 
see for example Horvath, Petrov, Scott and Showalter [TS] or Terman [29]. However 
our main interest here is the stability of these cellular fronts themselves. It has been 
shown by direct numerical simulation of the reaction-diffusion system that the cellular 
fronts become unstable as the autocatalyst diffusion parameter is increased further — 
see for example Horvath, Petrov, Scott and Showalter [18] and Malevanets, Careta 
and Kapral [23|. We apply our multi-dimensional Evans function shooting method 
to this problem and we reveal the explicit form of the spectrum associated with such 
fronts. We compare our shooting methods with standard direct projection methods 
obtaining high accuracy concurrency. We confirm the literature on these instabilities 
and establish the shooting method as an effective and accurate alternative to standard 
direct projection methods. 

Our paper is structured as follows. In Section [2] we review Grassmannian and 
Stiefel manifolds. We show how to construct the stable and unstable solution sub- 
spaces of the spectral problem using a given Grassmannian coordinatization and define 
the Evans function using coordinatizations for each subspace. We show how to project 
the spectral problem onto a finite dimensional transverse basis in detail in Section [3] 
A large portion of numerical effort is devoted to accurately constructing the multi- 
dimensional fronts whose stability we wish to study, and in Section |4] we detail the 
technique we used to do this. In Section [5] we outline our main application — the cubic 
autocatalysis problem — and review the known travelling wave solutions and stability 
properties. We then bring to bear our full range of numerical techniques to the cubic 
autocatalysis problem in Section [G] comparing their accuracy. In Section [7] we provide 
some concluding remarks and outline future directions of investigation. 

2. Review: Grassmann flows and spectral matching. 

2.1. Grassmannian coordinate representations. A m-frame is a m-tuple of 
m < n linearly independent vectors in C". The Stiefel manifold V(n, m) of m- frames 
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is the open subset of C"^™ of all m-frames centred at the origin. Hence any element 
Y G V(n, m) can be represented by an n x m matrix of rank m: 



(2.1) 



The set of m dimensional subspaces of C" forms a complex manifold Gr(7i, m) called 
the Grassmann manifold of m-planes in C" ; it is compact and connected (see Steen- 
rod [151 p. 35] or Griffiths and Harris [TH p. 193]). There is a quotient map from 
V(n, m) to Gr(n, m), sending each m-frame centred at the origin to the m-plane it 
spans — see Milnor-Stasheff [23J p. 56]. Any m-plane in C" can be represented by 
an n X m matrix of rank m, like Y above. However any two such matrices Y and 
Y' related by a rank m transformation u £ GL(A:), so that Y' — Yu, will represent 
the same m-plane. Hence we can cover the Grassmann manifold Gr(n, m) by coordi- 
nate patches Ui, labelled with multi-index i — {ii, . . . , im} C {1, . . . , n}, where Uj is 
the set of m-planes Y E Gr(n, m), which have a matrix representation j/io £ C"^™ 
whose m X m submatrix, designated by the rows i, is the identity matrix. Each such 
patch is an open dense subset of Gr(7i, m) isomorphic to For example, if 

i = {1, . . . , m}, then V{i....^rn} can be uniquely represented by a matrix of the form 



- [^;) (2.2) 

with 1° = {m + 1, . . . ,n} and y G c(n-ni)m each i, there is a bijective map 

ip\ : Ui^C*^""'"-''" given by ip\ : yjo i-^- y, representing the local coordinate chart for 
the coordinate patch Ui. For more details see Griffiths and Harris [Ml p. 193-4]. 

2.2. Riccati flow. Consider a non- autonomous linear vector field defined on the 
Stiefel manifold of the form V{x, Y) — A{x) Y, where Y e V(n, m) and A: R^gl(n), 
the general linear algebra of rank n matrices. Following the exposition in Ledoux, 
Malham and Thiimmler 21], fix a coordinate patch Ui for Gr(n, m) for some i, and 
decompose Y E V(n, m) into 

Y = j/iou. 

If we substitute this form into the ordinary differential system Y' = V{x,Y) we obtain: 

y[oU + yiou' = {Ai + Aioy^o)u, 

where A^ denotes the submatrix obtained by restricting the matrix A{x) to its ith 
columns. If we project this last system onto its i°th and ith rows, respectively, we 
generate the following pair of equations for the coordinate chart variables y = ip\o y^o 
and transformations u, respectively: 

y ^c + dy-ya- yby (2.3a) 

and 

u' = ia + by) u. (2.3b) 



Here a,b, c and d denote the i x i, i x i° , i° x i and i° x i° submatrices of A, respectively. 
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We observe that the flow on the Stiefel manifold Y{n, m) generated by the hnear 
vector field V{x^ Y) can be decomposed into an independent Riccati flow in the coor- 
dinate chart variables of a fixed coordinate patch Ui of Gr(n,rn), and a slaved linear 
flow of transformations in GL(to). Hence we can reconstruct the flow on the Stiefel 
manifold generated by the linear vector field V{x,Y), by simply solving the Riccati 
equation above, and then if required, subsequently integrating the linear equation for 
the transformations. In general, solutions to the Riccati equation (|2.3ap can blow up. 
However, this is simply the manifestation of a poorly chosen local representative coor- 
dinate patch. This can be gleaned from the fact that simultaneously the determinant 
of u becomes zero — afterall the linear Stiefel flow, with globally smooth coefficients, 
does not blow up. To resolve this, we simply change patch and keep a record of the 
determinant of the patch swapping transformation. Further details and strategies can 
be found in Ledoux, Malham and Thiimmler 21j. 

2.3. Drury— Oja flow. The continuous orthogonalization method of Humpherys 
and Zumbrun [19] can be thought of as the flow on the Grassmann manifold corre- 
sponding to the linear vector fleld V, with the coordinatization evolving according to 
a given unitary flow (see Ledoux, Malham and Thiimmler [21j). Indeed their method 
is specified by a Drury-Oja flow for the n x m orthogonal matrix Q, together with 
the flow for the determinant of the matrix R, in the Qi^-decomposition of Y: 

Q' = {In - QQ^)A{x)Q, (2.4a) 
(det R)' = Tr(Qt A(a;)g) (det R). (2.4b) 

Here denotes the Hermitian transpose of Q. In practice the determinant of R is 
exponentially rescaled — see Humpherys and Zumbrun |19j . We also did not find it 
necessary to apply any of the stabilization techniques suggested therein. 

2.4. Spectral problem. Consider the linear spectral problem on M with spec- 
tral parameter A in standard first order form with coefficient matrix A{x; A) € C"^": 

Y' = A{x;\)Y. (2.5) 

We assume there exists a subdomain C C containing the right-half complex plane, 
that does not intersect the essential spectrum. For A € 17, we know that there exists 
exponential dichotomies on and R+ with the same Morse index m in each case (see 
Sandstede [25]). Hence for A e il, let Y~{x;X) € V(n, m) denote the matrix whose 
columns are solutions to (|2.5p which span the unstable subspace of solutions decaying 
exponentially to zero as x— s- — oo, and let Y'^{x; A) € V(n, n — ra) denote the matrix 
whose columns are the solutions which span the stable subspace of solutions decaying 
exponentially to zero as x— *■ -I- oo. The values of spectral parameter A for which 
the columns of Y~ and columns of Y^ are linearly dependent on R are pure-point 
eigenvalues. 

2.5. Evans determinant. The Evans function D{\) is the measure of the linear 
dependence between the two basis sets Y~ and Y'^: 

i?(A) = e-/^TrA(C;A)dc r+(x;A)), (2.6) 

(see Alexander, Gardner and Jones [Tj or Sandstede [25] for more details). Note that 
for the decompositions Y^ = y\o u_ and Y^ = yi<j^u+, assuming det m_ 7^ and 
det M+ 0, we have 

det(y^ F+) = det (yt^ • det w_ • det m+. 
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Let Yq' (A) denote the nxm matrix whose columns are the m eigenvectors of A(— oo; A) 
corresponding to eigenvalues with a positive real part. Analogously let Yf^{X) denote 
the n X {n — m) matrix whose columns are the (n — m) eigenvectors of A{+oo; A) 
corresponding to eigenvalues with a negative real part. 

Suppose we fix a coordinate patch for Gr(n, m) labelled by i_ (which has car- 
dinality m). We integrate the corresponding Riccati differential problem (|2.3a|) for 
y_ — ipi_ o Uio^ from x = —oo to a; = x*. The initial data oo; A) is generated 
by postmultiplying the i°_ x {l,...,m} submatrix of Yq^{X) by the inverse of the 
i_ X {1, . . . , m} submatrix of yg"(A) (i.e. perform the simple decomposition into 
and U- for Y(^{X), and use the i!L x {1, . . . , m} block of j/i" as the initial data for y_). 
Similarly we fix a patch for Gr(n, n — m) labelled by (with cardinality (n — m)) and 
integrate the Riccati differential problem (|2.3ap for y+ = ipi^ o y^o from x = +oo to 
X = Xf. We perform the analogous decomposition on Y^{X) with index to generate 
the initial data y^{+oo; A). If the solutions to both Riccati problems do not blow up, 
so that detM_(a;; A) ^ for all x G (— oo,a;,] and detM+(a;; A) ^ for all x e [x»,oo), 
then we define the modified Evans function, which is analytic in A, by 

^(A;^;,) EEdet(yio (x,;A) yio(x,;A)), (2.7) 

where a;* € M is the matching point. In our practical implementation in Section O 
we took the patches labelled by i- = {1, . . . , m} and i+ = {m + 1, . . . , n} for the left 
and right problems, respectively. We used the generally accepted rule of thumb for 
Evans function calculations: to match at a point roughly centred at the front. In our 
example application we took a;, = 50, and did not observe singularities in y~ or y'^. 
However generally, varying may introduce singularities for y~ or y^ in their re- 
spective domains (— oo, a;*] and [x^,, oo). Indeed this was observed in several examples 
considered in Ledoux, Malham and Thiimmler [21j . There, a robust remedy is pro- 
vided, which on numerical evidence, allows matching anywhere in the computational 
domain (and which in future work we intend to apply to our context here) . 

If we use the Humpherys and Zumbrun continuous orthogonalization approach 
then we define the Evans function to be 

£'Hz(A;a;,) = det(g"(x,; A) Q+(x,; A)) • deti?-(x,;A) • det i?+(a;,; A). (2.8) 

This is analytic in A when we include the det terms. 

2.6. Angle between subspaces. We can also measure the angle between the 
unstable and stable linear subspaces V(n, to) and Y{n, n—m), see Bjorck and Golub [7]. 
When the columns of the matrices Y^{x;X) and Y^{x;X) which span V(n, m) and 
V(n, n — to), respectively, are nearly linearly dependent the angle will be small. It 
is defined as follows — without loss of generality we assume to > n/2. Let Q~ and 
Q+ be unitary bases for Y{n,m) and Y{n,n — to). They are specified by the QR- 
decompositions Y~ = Q~ R~ and Y^ = Q^R~^. The cosine of the smallest angle 8 
between V(n, m) and V(n, n — m) is given by the largest singular value of the matrix 
(Q )^Q^- However when 6 is small it is more accurate to compute sm9, which is 
given by the smallest singular value of 

Hence if the singular values of this matrix are (7^, i = 1, . . . , n — m, define 

0(A; a;,) = arcsin min ct^, 

i—l....,n — k 

where again, a;, is the matching point. Eigenvalues correspond to zeros of 0{X; x*). 
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3. Projection onto a finite transverse Fourier basis. Consider the parabolic 
nonlinear system (|l.ip posed on the cylindrical domain M x T. Recall that we assume 
Dirichlet boundary conditions in the infinite longitudinal direction a; G R and periodic 
boundary conditions in the transverse direction y (z T (here T is the one-torus with 
period/circumference L). Suppose we have a travelling wave solution U = Uc{x,y) 
of velocity c satisfying the boundary conditions. The stability of this solution is 
determined by the spectrum of the linear differential operator 

C = BA + cd^+I)F{Uc), (3.1) 

where DF is the Jacobian of the interaction term. Hence the eigenvalue problem is 

BAU + cd^U + I)F{Uc)U = XU. (3.2) 

We wish to define and compute the Evans function for this spectral problem. The idea 
is to perform a Fourier decomposition in the transverse direction. This transforms the 
problem to a one-dimensional problem, for which we can apply the Evans function 
shooting approach. Hence suppose that we can write 

oo 

U{x,y)= Uk{x)e''y/\ (3.3) 

k— — oc 

where Uk & arc the Fourier coefficients (with N denoting the number of compo- 
nents of U) and L — 2ttL — this implicitly restricts our attention to perturbations U 
with period L. Similarly, we expand DF{Uc) as 

oc 

BF{U,{x,y))^ J2 i^fe(x)e^'=«/^, (3.4) 

k— — oo 

where the Fourier coefficients Dk G C^^^. Substituting our Fourier expansions for 
U and DF([/c) in p.3p and (|3.4p in the eigenvalue problem (|3.2p yields 



k— — oo 

oo 



k—-~oo k— — oo 

oo oo oo 

iky / L 



i=z — oo k- 

Reordering the double sum 



oo ^ oo oo 

Jky/L 



1 — — OQ k — — oo k — — oo iy'—~oo 

yields that for all A: G Z, we must have 

oo 



V — — 00 



In practical implementation, we consider only the first K Fourier modes. Hence we 
replace the above infinite system of equations by following finite system of equations 
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which we now also express in first order form: 

d^tJk = Pfc, (3.5a) 

K 

d^Pk^\B-^Uk + {k/LfUk-cB'^Pk- B-^Dk-,U„ (3.5b) 



-K 



for k = —K, —K + 1, . . . , K . This is a large system of ordinary differential equations 
involving the spectral parameter A, so we can apply the standard Evans function 
approach to it. In matrix form this is 



dx 



^3 (a;; A) 



A, 



P 



(3.6) 



where U is the vector of C^-valued components U-k, U-k+i, ■ ■ ■ , Uk and P is the 
vector of C^-valued components P-k, P-k+i, • ■ ■ , Pk- The lower right block matrix 
is A4 = —cB~^ (g) I2K+1- The lower left block matrix is given by ^3 = A3(A) + ^3(2;) 
where if we set Ek{X) = \B~^ + {k/Lyiff then 



/E^k{\) O 
O E.K+i{\) 



and 



Di Do 



2K D2K-I 



o 



O E+K{\)j 



D-2K \ 
D-2K+1 

Do J 



A similar Galerkin approximation for nonplanar fronts can be found in Gesztesy, 
Latushkin and Zumbrun [T^l p. 28,37]. 

4. The freezing method. To study the long-time behaviour of the wrinkled 
front solution we use the method of freezing travelling waves described in Beyn and 
Thiimmler [5] and analyzed in Thiimmler [30] . This method transports the problem 
into a moving frame whose velocity is computed during the computation. For the 
travelling waves whose stability is our concern here, the freezing method has two 
advantages over direct simulation: Firstly, we can do long-time computations of time 
dependent (for example oscillating) solutions. This is because the method picks out 
the correct speed for the moving frame in which the wave is roughly stationary (in 
particular it does not drift out of the chosen domain) . Secondly, the freezing method 
can be used to obtain a suitably accurate guess for the waveform, which can then be 
used to initiate a Newton solver that computes the travelling wave as a stationary 
solution in the co-moving frame. For completeness we will briefly review the idea of 
the freezing method here. Consider a parabolic nonlinear system of form (II. ip in a 
stationary frame of reference: 



dtU = BAU + F{U). 



(4.1) 
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If we substitute the ansatz U{x, y, t) — V{x — ^{t), y, t) into this equation, where ^{t) 
is the unknown reference position for a wave travelUng in the longitudinal x-direction, 
we generate the partial differential equation for V and 7': 

dtV = BAV + Y{t)d^V + F(y). (4.2a) 

In order to compensate for the additional degree of freedom which has been introduced 
by the unknown position 7, we specify the following additional algebraic constraint 

0= / {d.,V{x,y,t)f {V{x,y,t) -V{x,y,t)) dxdy. (4.2b) 

jRxT 

This constraint selects a unique solution from the one-parameter group of possible 
solutions generated by longitudinal translation. It is a phase condition that normalizes 
the system with respect to a given template function — for example a suitable choice 
is to set V{x, y, t) to be the initial data V{x, y, 0). Hence we have generated a partial 
differential algebraic system (|4.2p . for the unknowns V{x,y,t) and C,{t) = 7'(t), with 
initial data V{x,y,0) and C(0) within the attractive set of the stable travelling wave 
and its velocity. 

A stationary solution (V, C) = {Vc,c) of the partial differential algebraic sys- 
tem (|4.2p represents a travelling wave solution Uc{x,y,t) = Vc{x — ct,y) with fixed 
speed c to the parabolic nonlinear problem (|4.ip . If the travelling wave Uc is stable 
then the stationary solution (Vc, c) is also stable and the time evolution of partial dif- 
ferential algebraic system (14. 2p will converge to {Vc,c). Thus the numerical method 
consists of solving the partial differential algebraic system (I4.2p until the solution 
(V, C) has reached the equilibrium {Vc, c). Alternatively one can also solve the partial 
differential algebraic system (|4.2p until some time t = T when the solution is deemed 
sufficiently close to equilibrium and then use a direct method (e.g. Newton's method) 
to solve the stationary equation 

= B Av + 0d^v + F{v), (4.3a) 

0= / {d,V{x,y)f{V{x,y)-vix,y))dxdy, (4.3b) 

JRxT 

for (w, 9) starting with initial values vq{x, y) = V(x, y, T) and 9a = C,{T). 

5. Review: autocatalytic system. As our model application wc consider the 
cubic autocatalysis system 

ut — Au + cdxU — uu^, (5.1a) 
Vt=5Av + cdxV + uv'^. (5.1b) 

Here m(x, y, t) is the concentration of the reactant and v{x, y, t) is the concentration of 
the autocatalyst, in the infinitely extended medium [x, y) g R x T. In the far field, we 
suppose (u,w) approaches the stable homogencnous steady state (0, 1) as x ^ — cxi, 
and the unstable homogeneous steady state (1,0) as x — > -l-oo. The parameter 5 is 
the ratio of the diffusivity of the autocatalyst to that of the reactant. This system is 
is globally well-posed for smooth initial data, and any finite 6 > 0. 

From Billingham and Needham T we know that, in one spatial dimension, a 
unique heteroclinic connection between the unstable and stable homogeneous steady 
states exists for wavespeeds c > Cmin. The unique travelling wave for c = c,nin con- 
verges exponentially to the homogeneous steady states. They are readily constructed 
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by shooting (see for example Balmforth, Craster and Malham [3i) and represent the 
planar front solution referred to hereafter. 

It is well known that these planar travelling wave solutions become unstable 
to transverse perturbations when S is greater than a critical ratio Jcr ~ 2.3 (from 
Malevanets, Careta and Kapral [23l p. 4733]), and our transverse domain is sufficiently 
large. This is known as the cellular instability — see for example Terman [29], Horvath, 
Petrov, Scott and Showalter |l8j and Malevanets, Careta and Kapral [23]. Calculating 
the spectrum that reveals the mechanism of the instability is also straightforward. If 
we consider small perturbations about the underlying planar travelling wave solution 
of the form C/(x)e'''*+''^'"'/^)^, for any fc e Z, then the associated linear operator C 
whose spectrum we wish to determine has the form 

C = B{d^, - {k/Lfl) +cd^+ r)F{Uc). (5.2) 

The spectrum of £ consists of the pure-point spectrum, i.e. isolated eigenvalues of 
finite multiplicity, together with the essential spectrum. By a result of Henry [15j . if 
we consider £ as an operator on the space of square integrable functions, we know 
that the essential spectrum is contained within the parabolic curves of the continuous 
spectrum in left-half complex plane. The continuous spectrum for this model problem 
is given for all G R and each fixed fc € Z, by the locus in the complex A-plane of the 
curves 

A = ~5{k/LY + ic/i - 

A = -1 - 5{klLf + lc^l - 5fi^, 

A = -{k/iY +ic/i- /x^, 

where the last curve appears with multiplicity two. Hence the continuous spectrum 
touches the imaginary axis at the origin when k is zero. The origin is a simple 
eigenvalue for C due to translational invariance of the underlying travelling wave 
Uc{x) in the longitudinal direction. Thus, even though the essential spectrum does 
not affect linear stability, we cannot automatically deduce asymptotic stability unless 
we consider a gently weighted space of square integrable functions, which will shift 
the essential spectrum to the left. 

We now notice that the spectral problem for the linear differential operator £ is 
fourth order and one-dimensional, the transverse perturbation information is encoded 
though the transverse wavenumber parameter k. Hence we can employ the Evans 
function shooting approach, in particular in Figure [5T] we plot the dispersion relation 
for planar fronts corresponding to different values of the physical parameter 5. The 
dispersion relation gives the relationship between the real part of the exponential 
growth factor of small perturbations, Re(A), and transverse wave number k (in fact 
to make comparisons later easier we have used k/L on the abscissa). To construct 
the graph shown, we locate zeros of the Evans function on the real A axis and follow 
them as the transverse wave number k is varied. We see in Figure ISTTl that the interval 
(0, ks) of unstable transverse wavenumbers k increases as 5 increases (at least up to 
5 = 5). The wavenumber where the maximum of Re(A) occurs has largest growth and 
is thus typically the transverse mode expressed in the dynamics. 

However our real concern here is the stability of the wrinkled cellular fronts them- 
selves as we increase the diffusion ratio parameter 5 beyond 5c-c- We know the planar 
waves are unstable in this regime and the question is whether the bifurcating wrinkled 
cellular fronts are stable for values of 5 just beyond 5cr but do become unstable at 
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wave number /(/ (27t/.) 

Fig. 5.1. Dispersion relations for the planar front for different values of 5. For 5 = 2, all 
perturbations decay and the planar front is linearly stable, but for the other values of 5 there is a 
growing interval (0, k^) of unstable transverse wavenumbers. 

another critical value of 6. Using direct simulations of the reaction-diffusion system 
Malevanets, Careta and Kapral [53] demonstrated that when 6 = 5 the transverse 
structure of the propagating fronts was very complex with spatio-temporal dynamics 
of its own. Our goal is to elicit in practice the initial mechanisms that precipitate 
the behaviour they observe through a careful study of the spectra of the underlying 
multi-dimensional travelling fronts. 

6. Numerics. 

6.1. The method. We will compute the spectrum of the two-dimensional trav- 
elling front by computing the Evans function associated with the linear operator C 
defined in (|3.ip . First, we compute the front with the freezing method as explained 
in Section [4l Then as in Section [3l we project the problem onto a finite number of 
transverse Fourier modes k = —K, . . . , K . This generates the large linear system (j3.6p 
which is of the standard linear form (|2.5p where the coefficient matrix is given by 

•^(-^^) = fe," '"T'')- 

where the matrices and A4 are explicitly given at the end of Section [S] This 
coefficient matrix depends on both the travelling wave solution Uc and the spectral 
parameter A. It also depends on the nonlinear reaction term. For our autocatalytic 
problem the Jacobian D_F of the nonlinear reaction term and its Fourier transform D 
are given by 

BF{u) = (-f £)=(-!'*!' 

^ ' y V 2uv J \^ V * V 2u* V J 

where u and v are the Fourier transforms of u and v respectively, and * denotes the 
convolution defined by (u * v)k = J2tL-oo ^f^k-i- 

As described in Section [5751 we construct the unstable subspace Y^^ {X) of p.6p at 
X = —00, as the matrix whose columns are the eigenvectors of A{—oo; A) corresponding 
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to eigenvalues with a positive real part. Similarly we construct the stable subspace 
Yf^{X) at X = +00 from the eigenvectors of yl(+oo; A) corresponding to eigenvalues 
with a negative real part. Then to construct the modified Evans function (12. 7p at the 
matching point we proceeded as follows (for both the planar and wrinkled front 
cases for the autocatalytic problem below). The full system is of size n ~ 4(2X + 1). 
We use the coordinate patch identified by i_ = {1, . . . ,m} where m — 2{2K + 1) 
for the interval (— oo,a;*]. On this interval we solve the Riccati equation ()2.3ap for 
with the chosen coordinate patch, the coefficient matrices in the Riccati equation 
are a — Om, b — Im, c — A3 and d — A^. The initial data for j/_ is simply the 
lower (n — m) x m block of (A) postmultiplicd by the inverse of the upper m x m 
block of yQ^(A). For the interval [x^, +00) we use the coordinate patch identified by 
i+ = {m + l,...,n}. We solve the Riccati equation (|2.3ap for the coefficient 
matrices in this case are a = A4, b = A^, c = Im and d — 0,n- The initial data for 
is the upper m x (n — m) block of Yf^{X) postmultiplied by the inverse of the lower 
(n — m) X (n — m) block of Fq^(A). Of course in practical calculations, our intervals of 
integration are and [a:*,L^] for suitable choices of L~ and L+. In particular, 

we assume the interval L+] contains the travelling front and extends far enough 
in both directions that we can assume the front is sufficiently close to the far field 
homogeneous steady states. Further, the respective coordinate patches for the left and 
right intervals, and matching point chosen to avoid singularities in the Riccati 

flows in these intervals (in all our calculations for the autocatalysis problem we were 
always able to find such a matching point when using these two patches). Hence we 
evaluate the Evans function by computing the determinant (12. 7p . This information is 
used to find the spectrum of the travelling wave. 

6.2. The planar front. We first consider planar travelling waves. This will not 
give us new information but it is a good test for our procedure. The travelling waves 
can be found by the freezing method, but also by the shooting method explained 
in Balmforth, Craster and Malham [3] for the one-dimensional case. The next step 
is to find the unstable subspace of ()3.5p at x — —00. The travelling wave Uc is 
planar and thus independent of y, so the Fourier coefficient Dk vanishes for k ^ Q. 
This implies that the differential equation (|3.5p decouples in 2K + 1 systems, one for 
every k = —K, . . . ,K, each of the form 



We have (mi,M2) (0, 1) in the limit x — > —00. The coefficient matrix in (|6.ip has 
two unstable eigenvectors when A is to the right of the continuous spectrum, and these 
eigenvectors are 



9xUk — Pk, 



For the autocatalytic system (15. ip . this becomes 

/ 1 ° \ 

^ (^^\ \ 1 




(6.1) 




(6.2a) 
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and 

{0,1,0, fif, where = -^c + ^ + (i-f + X. (6.2b) 

These eigenvectors are ah analytic functions of A except at those points where either 
v or one of the expressions under the square root sign vanishes. A smaU calculation 
shows that the latter can only happen when A is on the negative real axis, so it will 
not interfere with stability computations, which concern eigenvalues with positive real 
part. In contrast, we may have u — for positive A, but in the situations considered 
here v vanishes only for values of A which are so large that they do not concern us 
here. Alternatively, the singularities caused hy v — can easily be avoided by using 
{v, 1, fjiv, /i)""" as eigenvector instead of (1, l/v, fi, ^i/v)^ . 

We thus have two unstable directions for every wavenumber k. Together these 
form the 2{2K + l)-dimensional unstable subspace Yq {\) at a; = — oo. Using this, as 
described in Section [Ql above, we integrate the Riccati equation (|2.3a|l from x = L~ 
to X = x^, where corresponds to a position sufficiently far behind the front. The 
integration is done with Matlab function ode45, which is based on the explicit fifth- 
order Runge-Kutta method due to Dormand and Prince. In all experiments, we use 
an absolute tolerance of 10"^ and a relative tolerance of 10~^. 

In the other limit, where x — s- -foo, we need the two stable eigenvectors of the 
coefficient matrix in (j6.ip . which are 

(1,0,M,0)T, where A' - - V^i^ + (|)' + i (6.3a) 

and 

{0,1,0, where ^^ ^ - ^J {d^ + {ff + X. (6.3b) 

Again, we get a 2{2K + l)-dimensional subspace {X). We use this to generate the 
initial condition for the Riccati equation (|2.3ap we integrate backwards from x = to 
X = — note the coefficients are different to those for the left interval — see Section l6T] 
above. Here is a longitudinal position sufficiently far ahead of the front. In our 
actual calculations we used — ±25 and a;* = 0. Finally, the Evans function is 
computed using the modified form (j2.7[) . 

This concludes the explanation for the Riccati approach. The process for the 
Drury-Oja approach due to Humpherys and Zumbrun T9' is very similar. The dif- 
ference is that it uses the Drury-Oja flow (|2.4p to propagate the subspaces and (|2.8p 
to evaluate the Evans function. 

As an illustration, we show in Figure 16.11 a graph of the Evans function along the 
real axis for (5 = 3, with the Fourier cut-off dX K = 24. We see that the Evans function 
has a zero at A = and double zeros around 0.0005, 0.0007 and 0.0011. The zero at 
A = corresponds to the eigenvalue at the origin caused by translational symmetry. 
The other zeros correspond to double eigenvalues. They have positive real part and 
thus we can conclude that the planar wave is unstable for 5 = 3 — coinciding with the 
long-established statements in Section [5l 

In fact, the zeros for the Evans function can be related to the dispersion curve in 
Figure [5TTI As mentioned above, the differential equation ()3.5|) decouples into 2K -\- 1 
subsystems of the form (|6.ip , each corresponding to the linear operator C given in (|5.2p 
for a particular wave number. The basis vectors we chose for the unstable subspace 
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Fig. 6.1. The Evans function along the real axis for the planar front at 5 = 3. The Fourier 
cut-off is K = 2A and the transversal length of the domain is L = 200. The Evans function is 
computed using the Drury-Oja approach. A very similar plot can be produced using the Riccati 
approach. 

are zero for all but one of these subsystems. Since the subsystems are decoupled, the 
solution bX X — is also zero for all but one of these subsystems. Thus, the matrix 
{Y~ {x*]X) y"'"(a;*;A)) at the matching point is (after reordering) a block-diagonal 
matrix consisting of 2K + 1 four-by-four blocks, and its determinant is the product 
of the determinants of the 4x4 blocks. However, every determinant of a sub-block 
is the Evans function of the operator (|5.2p for a particular wave number k. Thus, 
the Evans function for the two-dimensional planar wave is the product of the Evans 
functions for the one-dimensional wave with respect to transverse perturbations: 

K 

D{\)^ n ^''W- (6-4) 

k=-K 

We chose L — 200, so fc = 1 corresponds to a wave number of ^ w 0.0314 in Fig- 
ure [O] The plot shows that the corresponding growth rate is approximately 0.0005, 
thus Di{X) has a zero around A — 0.0005. The dispersion relation is symmetric, so 
Z)_i(A) has a zero at the same value. This explains why the two-dimensional Evans 
function has a double zero around A = 0.0005. Similarly, the double zeros around 
A = 0.0011 and A = 0.0007 correspond to fc = 2 and fc = 3, respectively. 

As Figure 16.11 shows, the Evans function is of the order 10~^° in the interval of 
interest. It is important to emphasize that we are primarily interested in the zeros 
of the Evans function, which correspond to eigenvalues of the operator C defined 
in p.ip . The scale of the Evans function between the zeros, which depends on the 
normalization of the eigenvectors (|6.2p and (|6.3p . is of lesser relevance in this context. 
However for completeness, we address this issue in Section [6.51 

6.3. The wrinkled front. Having established that our method works for planar 
fronts, we now turn to the wrinkled front. The most immediate problem is that 
of computing reliable steady wrinkled fronts. The procedure we used to do this is 
outlined as follows: 
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Fig. 6.2. The wrinkled front for 5 = 3, The left panel shows the u component and the right 
panel shows the v component. 

• The partial differential algebraic system (|4.2p is solved on the computational 
domain [—150,150] x [—60,60] with a second-order finite element method 
on a grid of right triangles formed by dividing a grid of 300 x 240 rect- 
angles along their diagonals. We used the Finite Element package Comsol 
Multiphysics"^'^ [8] which includes a version of the DAE solver daspk. 

• To obtain reasonable starting waveform profiles for the partial differential 
algebraic system (j4.2p , we constructed the one-dimensional waveform, which 
we then swept out uniformly in the transverse direction. 

• For travelling waves that appeared to be stable, we use the freezing method 
described in Section 2] to obtain good starting waveform profiles to initiate 
the use of Newton's method in Comsol. 

• For travelling waves that appeared to be unstable, we use simple parame- 
ter continuation (also using Comsol) starting with a (parametrically) nearby 
stable wave. 

Figure shows the front which results for (5 = 3. 

The wrinkled front varies along the transverse y-direction. Therefore, the Fourier 
coefficients Dk do not vanish for fc 7^ and the differential equation p.6p does not 
decouple. In the limits x — s- ±00, we have restricted ourselves to wrinkled front profiles 
that approach a y-independent state, which is the same as for the planar front. Hence, 
the computation of the stable and unstable subspaces for the planar front, which we 
reported above, remains valid for the wrinkled front. 

The computation now proceeds as with the planar front. We solve the Riccati 
differential equation (|2.3ap in the left and right intervals (as described in Section WA] 
above) for the Riccati approach, and (|2.4p for the Drury-Oja approach. The front Uc 
is determined as explained above, with interpolation being used for those points that 
fall between the grid points of the finite clement solution. We use cubic interpolation, 
because linear interpolation is not sufficiently accurate. Finally, the Evans function is 
computed using either the modified form (j2.7p or Humpherys and Zumbrun form ()2.8p . 

Figure 16.31 shows a plot of the Evans function along the real axis for 6 = 3, 
computed using the Riccati approach. The result of the Drury-Oja approach is shown 
in Figure 16.41 For both approaches, we used projection on K = 9 modes and we 
solved the Riccati differential equations ()2.3ap (with appropriate coefficients for the 
left and right intervals), and the Drury-Oja equations (|2.4I) . by the Matlab ODE- 
solver ode45 with absolute and relative tolerances of 10^^ and 10~^. The domain is 




Fig. 6.4. Same as Figure R?. 3\ but with the Drury-Oja approach. 
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Eigenvalues (Evans 
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3 


0.001609 


-0.000026 


-0.000781 


-0.001296 


-0.000670 


-0.006078 


-0.006177 


4 


0.001609 


0.000002 


-0.000001 


-0.000519 


-0.000670 


-0.006078 


-0.006177 


rr 




0.001589 


0.000002 


-0.000001 


-0.000519 


-0.000720 


-0.005295 


-0.005426 


D 


0.001589 


-0.000002 


-0.000003 


-0.000515 


-0.000720 


-0.005295 


-0.005426 


7 


0.001589 


-0.000002 


-0.000003 


-0.000515 


-0.000721 


-0.005292 


-0.005422 


8 


0.001589 


-0.000002 


-0.000003 


-0.000515 


-0.000721 


-0.005292 


-0.005422 


9 


0.001589 


-0.000002 


-0.000003 


-0.000515 


-0.000721 


-0.005292 


-0.005422 


24 


0.001589 


-0.000002 


-0.000003 


-0.000515 


-0.000721 


-0.005292 


-0.005422 




Eigenvalues (ARPACK) 




0.001592 


0.000000 


0.000000 


-0.000514 


-0.000719 


-0.005289 


-0.005420 



Table 6.1 

Zeros of the Evans function for the wrinkled front at 5 = 3 for different values of the wave 
number cut-off K , as computed with the Riccati and Drury-Oja approaches. The last row shows the 
eigenvalues computed by ARPACK. 



[—150, 150] X [—60, 60] and as matching point we have chosen a;, — 50 which is roughly 
the front location (see Figure 16. 2p . Both approaches produce Evans functions with 
zeros at (approximately) the same A-values: there is one eigenvalue with a positive 
real part around A = 0.0016, a double zero around (due to translation invariance), 
and pairs of eigenvalues around —0.0006 and —0.0053. In particular, the eigenvalue 
around A = 0.0016 shows that the wrinkled front is unstable for 6 — 3. 

However, the Riccati and Drury-Oja approaches do differ in some respects. The 
Riccati approach involves integrating a system of half the dimension as compared to 
the Drury-Oja approach, so we can expect it to be faster, especially for large values of 
K (in comparable experiments for large K, it was typically two to three times faster). 
See Ledoux, Malham and Thiimmler |21j where some detailed accuracy versus cost 
comparisons can be found. We also note that the Drury-Oja approach again produced 
very small values for the Evans function, which as explained in Section is not as 
alarming as it first seems. The Riccati approach, which uses a different scaling of the 
Evans function, yields only moderately small values, though again see Section [6.51 

We also used the angle between the unstable and stable subspace 0(A;a;*) to 
determine the location of the eigenvalues. Though this scales better than the Evans 
determinant, it is non- negative by definition and zeros occur as touch-downs in the 
complex spectral parameter plane. Actual zeros of the function are thus harder to 
establish compared to sign changes in the real and imaginary parts of the Evans 
determinant. Further, 9{X;x^) is not analytic in A. 

More accurate values for the zeros of the Evans function can be found by using 
a standard root-finding method. This yields the values listed in Table 16.11 The 
results of the Riccati and Drury-Oja approaches agree (to the precision shown) in 
all cases. The table also shows results for different values of the wave number cut- 
off K . The eigenvalues converge quickly when K increases, and projection on K = 7 
modes is sufficient to get six digits of accuracy, at least in this example. The largest 
computation we performed is with K = 24, in which case the one-dimensional equation 
has order 196. 
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Fig. 6.5. The locus where the real and imaginary parts of Evans function vanish for the virinkled 
front at 5 = 3. This plot uses a 100 X 100 grid in the depicted region of the complex plane. The 
Evans function was computed on the 5000 grid points in the top half and extended by symmetry to 
the bottom half. 



We also computed the eigenvalues by a projection method to check these results. 
The Comsol package, which we used to compute the travelling front, has an eigenvalue 
solver based on the Arnoldi method as implemented in ARPACK [22]. This solver 
uses a shift a to allow the user to choose on which part of the spectrum to concentrate. 
We found that the choice of this shift is not straightforward. It is unwise to use a 
shift which is exactly an eigenvalue, so in particular, the default value of u = should 
not be used. We used a shift of cr = 1 or cr = 3 which proved to be safe. The 
resulting eigenvalues are also listed in Table 16.11 They are in close agreement with 
the eigenvalues found using the Evans function, giving additional confidence in our 
computation. 

Figure Wl5\ shows the contours where the real and imaginary parts are zero; eigen- 
values correspond with the intersections of these contours. All eigenvalues are real. 
Interestingly, the contours are well separated on the boundary of the depicted re- 
gion of the complex plane even though the eigenvalues are packed closely together. 
This suggests that it may be better to study the Evans function in a contour sur- 
rounding the eigenvalues instead of concentrating on the eigenvalues themselves — see 
Section EJl 

Figures 16.61 and 16.71 show the Evans function for the case S — 2.5. The Evans 
function no longer has a zero with positive real part, meaning that there are no 
unstable eigenvalues. Hence for this value of S the wrinkled front is in fact stable. 
We have thus established the following overall stability scenario for travelling wave 
solutions of the cubic autocatalysis problem. When S is less than da ~ 2.3 planar 
travelling waves are stable to transverse perturbations. For S beyond this instability 
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Fig. 6.6. The Evans function along the real axis for the wrinkled front at 5 = 2.5, computed 
using the Riccati approach. The right panel zooms in on the fiat plateau around the origin in the 
left panel. 




X 10 X 10 

Fig. 6.7. Same as Figure Y6. 61 but with the Drury-Oja approach. 



there exist wrinkled fronts that are stable for S — 2.5. However for (5 = 3 these 
wrinkled fronts themselves become unstable. These stability results may depend on 
the transverse domain size (which we fixed at L = 120); see Horvath, Petrov, Scott 
and Showalter [18]. 

6.4. Global eigenvalue search methods. If one is only interested in the sta- 
bility of the underlying front then several methods exist to detect eigenvalues in the 
right-half complex spectral parameter plane — which also require relatively few Evans 
function evaluations. One important method uses the argument principle. We have 
the following sectorial estimate for the spectrum of £: if A G cr(£) then 

Re(A) < ||DF(C/e)||Lo=(RxT), (6.5a) 

Re(A) + |Im(A)| < £^ + 2 ||DF(C/,)||i^(RxT), (6.5b) 

where k = min{i3ii, B22} — we provide a proof in Appendix \X[ Hence we can restrict 
our search for unstable eigenvalues to the region of the complex A-plane bounded by 
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Fig. 6.8. The left figure shows the contour C in the complex plane. The right figure shows 
the argument of D{X) when X transverses the top half of this contour (for the wrinkled front with 
5 = 3). The plot is split up in four parts, corresponding to the quarter circle, the segment along 
the imaginary axis, the diagonal segment, and the segment going down. The scale on the horizontal 
axis is linear on the first, third and fourth part and logarithmic on the second part. 

these inequalities. 

In our case, we have c ~ 0.577 and ||DF(C/c)||l°' ~ 1.422 for S = 2.5, and 
c « 0.548 and |jD_F(J7c)||L°' ~ 1.450 for 6 = 3. So, the above estimate imphes 
that any eigenvalue A satisfies Re(A) < 1.5 and Re(A) + |Im(A)| < 3. Any unstable 
eigenvalue is thus inside the contour bounded by these lines and the imaginary axis; 
this contour is depicted in the left part of Figure 16.81 (the small semi-circle around 
the origin is to exclude the double eigenvalue at the origin, as explained below). The 
Evans function is analytic inside this contour, and thus we can count the number 
of zeros inside the contour by computing the change of the argument of the Evans 
function as A goes around the contour. 

This suggests the following method. We divide the contour C in small segments, 
compute the Evans functions at the end points of the segments, find the change in 
the argument over each section and sum these to find the total change over the whole 
contour. The number of zeros is then the change of the argument divided by 27r. The 
segments in which the contour is divided have to be chosen small enough that the 
change of the argument over each segment is smaller than tt. Adaptive procedures 
to achieve this have been proposed by, amongst others, Ying and Katz [31], but for 
simplicity we determined the partition of the contour by trial and error. 

The Evans function, with our choice of initial conditions, satisfies D{X) = D{X). 
Hence, it suffices to transverse half the contour C. The plot at the right of Figure 
shows how the argument of D{X) changes as A transverses the top half of the contour. 
We see that the change over the whole contour is —2tt, and thus there is one zero 
inside the contour. This is the zero around A = 0.0016 shown in Figure there are 
no other unstable eigenvalues. 

The operator C is invariant under translations and hence it has two eigenvalues 
at A = 0. Thus, the Evans function is zero at A = and its argument is undefined. 
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For this reason, the contour C includes a smah semi-circle of radius 10""* to avoid 
the origin. It is of course possible that we miss some eigenvalues because of this 
deformation. To guard against this, we also compute the change in the argument 
of D{X) as A transverses a circle around the origin with radius 10"**. We find that 
the argument changes by 47r and hence the Evans function has two zeros inside the 
circle. This shows that we have not missed any eigenvalues by the deformation of C 
around the origin. 

Another method is the parity-index argument (which we have not investigated 
here, but we include a description for completeness). To determine if there are eigen- 
values on the positive real axis, we could use that 

^(A) - 26-^'-^/^, 

as A ^ +00, as we prove in Appendix |B] A change in sign of the Evans function along 
the real axis indicates an unstable eigenvalue. This can be exposed by comparing the 
sign of the Evans function in the asymptotic limit A — s- -l-oo to its sign elsewhere on 
the positive real axis, or the derivative or second derivative of the Evans function 
along the real axis evaluated at the origin. Note that to implement this technique 
we need to ensure the Evans function evaluated in the neighbourhood of the origin 
and that generating the asymptotic limit have been normalized, as far as their sign is 
concerned, in the same way. This may require keeping careful track of the the non-zero 
analytic multiples dropped in our definition of the modified Evans function, i.e. the 
normalization of the initial data, the scalar exponential factor and the determinants 
of u_ and u+ (which can be computed analytically in the asymptotic limit also) . 

6.5. Convergence of the Evans function. The results reported in Table [iT] 
show that the values of A for which the Evans function D{X) is zero converge as 
K ^ oo. However they do not show how the Evans function scales for other values of 
A (not in the spectrum) as K becomes large. We see in Figure I5TT] for the planar front, 
that in the domain of interest with A of order 10"'*, the Evans function computed 
using the Drury-Oja approach with K = 24, is of order 10~^°. Similarly, for the 
wrinkled front, we see in Figures 16.41 and 16.71 also computed using the Drury-Oja 
approach with K = 9, the scale of the Evans function is of order 10~^^ for A of order 
10~^. These small results we find for the Evans function — far smaller than machine 
precision — suggest that the computation might be inaccurate because of round-off 
error. This is not the case here. The solution of the Drury-Oja equation (|2.4[) does 
not contain such small numbers. They only appear when the determinant is computed 
at the matching point, which is a stable computation; see for example Higham [T71 
§13.5] . The determinant is computed by performing an LU-decomposition with partial 
pivoting and then multiplying the diagonal entries. Both steps are numerically stable 
even though the final result may be very small (as long as it does not underflow). 
The factorization (j6.4p for the planar front provides another view point: the Evans 
function for if = 24 is the product of 2K -I- 1 = 49 one-dimensional Evans functions, 
and if each of them is 0.06 — not a very small number — then their product is of the 
order 10"^", which is a very small number. 

For the modified Evans function wc construct using the Riccati approach, we see 
in Figure 16.31 that in the domain of interest, it is of order 10~^ when K = 9. If 
we reproduce Figure [6?3l for increasing values of K, we find that the modified Evans 
function actually grows. For large K, it increases roughly by a factor of 10 per unit 
increase in K. However it is important to emphasize that the solutions to the Riccati 
systems, that are used to construct the modified Evans ftmction, remain bounded. It 
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is only when we evaluate the determinant, that for example when K — 21, the scale 
of the modified Evans function in Figure [^751 is of order 10^. Hence the question of 
the scaling of the Evans function comes into play in the multi-dimensional context. 

In the one-dimensional case, the standard definition of the Evans function is (|2.6p 
given in Alexander, Gardner and Jones [1]. Of course this definition stands modulo 
any analytic non-zero multiplicative factor — afterall it's the zeros of the Evans func- 
tion we are typically after. A different appropriate normalization of the initial data 
yg^(A) rescales the Evans function by such a factor. Further it is common to drop the 
scalar exponential factor. Hence for example the difference between the Humpherys 
and Zumbrun Evans function (|2.8p and the standard Evans function is that the scalar 
exponential factor is dropped and the initial data might be normalized differently 
when it is Qi^-factorized. For the modified Evans function in (j2.7p . we dropped the 
exponential factor and also the determinants of some full rank submatrices in the 
determinant of the standard Evans function. 

However such analytic non-zero factors impact the scale of the Evans function 
in the multi-dimensional context, as they depend on K. This can lead to growth or 
decay of the Evans function employed away from the eigenvalues as K is increased, 
depending on which factors are kept and which are dropped. This is precisely what 
we have observed for the Drury-Oja and Riccati approaches, which decay and grow, 
respectively, as K increases. 

Gesztesy, Latushkin and Zumbrun [TU §4.1.4] show that the Evans function, when 
suitably normalized, coincides with the 2-modified Fredholm determinant of a certain 
operator related to C. Importantly, the 2-modified Fredholm determinant converges as 
K oo. They treat the self-adjoint case; the non-self-adjoint operator C we consider 
can, using suitable exponential weight functions, be transformed into a self-adjoint 
one (numerically though this is not necessary and could introduce further complica- 
tions due to the introduction of exponential factors in the potential term). Gesztesy, 
Latushkin and Zumbrun define the Evans function much like we have, through a 
Galerkin approximation, using the matching condition (j2.6l) . Their normalization, 
natural to their context, involves ensuring the Evans function does not depend on the 
coordinate system chosen and is invariant to similarity transformations (see Gesztesy, 
Latushkin and Makarov [llj). Crucially, Gesztesy, Latushkin and Zumbrun demon- 
strate in Theorem 4.15 that their 2-modified Fredholm determinant equals the Evans 
function multiplied by an exponential factor (nonzero and analytic), whose exponent 
diverges as K oo (it can diverge to positive or negative infinity depending on the 
precise form of the coefficient matrix A{x; A)). Hence their Evans function diverges 
as well. Consequently, the 2-modified Fredholm determinant is a natural candidate 
for the canonical Evans function. 

With this in mind, our observation regarding the measured divergence of the 
modified and Humpherys and Zumbrun Evans functions in our experiments, is not 
too surprising. Their divergence, the modified Evans function to large values, and 
the Humpherys and Zumbrun Evans function to small values, is simply a reflection of 
the over/under estimate of the crucial, correct multiplicative exponential factor that 
equates these Evans functions with the 2-modified Fredholm determinant of Gesztesy, 
Latushkin and Zumbrun, which guarantees a finite limit a,s K oo. It would be 
interesting to see whether one can explicitly compute the factor relating the Evans 
function of Gesztesy, Latushkin and Makarov and the Evans function, as defined in 
this paper. 
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7. Concluding remarks. Our goal in this paper was to show, for the first time, 
that spectral shooting methods can be extended to study the stability of genuinely 
multi-dimensional travelling fronts. It was already obvious to the people working in 
this field that this can be done by a Fourier decomposition in the transverse dimen- 
sion, thus reducing the problem to a large one-dimensional problem; however, no-one 
had yet managed to tackle the resulting large one-dimensional problems. This issue 
was overcome in theory by the work of Humpherys and Zumbrun |19j and latterly 
that of Ledoux, Malham and Thiimmler [5T] , which showed that enough information 
for matching can be retained by integrating along the natural underlying Grassman- 
nian manifolds (though these ideas have been well known in the control theory and 
quantum chemistry for a long while — see for example Hermann and Martin [16j and 
Johnson [20J). The dimension of the Grassmannian manifolds scale with the square 
of the order of the original problem. These methods brought the possibility of multi- 
dimensional shooting within our grasp, but there were still challenges to be overcome. 
We have addressed the most immediate of these in this paper: 

1. How to project the problem onto the finite transverse Fourier basis and ensure 
a numerically well-posed large system of linear ordinary differential spectral 
equations; 

2. How to construct the genuinely multi-dimensional fronts, this was a particu- 
lar challenge which we managed to overcome using the combined techniques 
of freezing and parameter continuation (the latter for the unstable multi- 
dimensional fronts); 

3. Choose a non-trivial but well-known example problem, which exhibits the de- 
velopment of planar front instabilities (an important initial testing ground for 
us) but also two-dimensional wrinkled fronts that themselves became unstable 
as a physical parameter was varied; 

4. Apply the Humpherys and Zumbrun continuous orthogonalization as well as 
the Riccati methods to study the stability of the wrinkled fronts, in particular 
integrating a large system of order 196 — using K — 24 transverse modes — 
hitherto the largest system considered in Evans function calculations was of 
order 6 (see Allen and Bridges [2]). 

5. Ultimately demonstrating the feasibility, accuracy, numerical convergence 
and robustness of multi-dimensional shooting to determine eigenvalues. 

However there is still more to do. Though multi-dimensional shooting is competitive 
in terms of accuracy, how does it compare in terms of computational cost? There are 
multiple sources of error including: approximation of the underlying wrinkled front, 
domain truncation, transverse Fourier subspace finite projection, approximation of 
the solution to the projected one-dimensional spectral equations and the matching 
position. Our experience is that the accuracy of the underlying wrinkled front is 
especially important. Although we have endeavoured to keep all these errors small, not 
only for additional confidence in our numerical results, but also for optimal efficiency, 
we would like to know the relative influence of each one of these sources of error in 
order to target our computational effort more effectively. Finally we would also like 
to see how the method transfers to hyperbolic travelling waves and more complicated 
mixed partial differential systems. 
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Appendix A. Sectorial estimate on the spectrum. We follow the standard 
arguments that prove that the parabolic linear operator £ is sectorial, see for example 
Brin [F, p. 26]. First consider the complex inner product with the spectral problem 

BAU + cd^U + BF{Uc) U = XU. 

Premultiplying this spectral equation with [/^, and integrating over the domain R x T, 
we get 

-[ \/U^B\/Udxdy + c I d^U + U'' J:)F{Uc)U dxdy ^ X\\U\\l2. 

In this last equation, if we set B = P^P where P = diag{i?J(^, i?22^}i S^^ 

MlUWh + WP'^UW^^c [ U'^ d^U + U^V)F{Uc)Udxdy. (A.l) 

JKxT 

If we consider the complex conjugate transpose of the spectral equation, postmultiply 
with J7, and integrate over the domain M x T, we get 

A||C/||^2 + ||PVC/||i2 = c / d^U^ U + {BFiUcif U dxdy. 

If we now add these last two equations then we get 

MmU\\h + \\PVU\\l,^ f U^^{DF+{DFf){U,)Udxdy 

< |!Di^(C/,)|U^||;7||i.. (A.2) 

Now dividing through by ||C/||^2 we get (|6.5ap 

Next consider taking the imaginary part of (|A.1[) followed by the absolute value 
on both sides and using Holder's inequality to get 

|Ini(A)| < c||C/|U2||VC/|U2 + ||DF(C/,)||i^||[/||2,. (A.3) 

Adding inequalities (IA.2P and (jA.SP and using that ||PV?7||^2 > K||Vt/||^2, where 
K = minjBii, B22}, and the arithmetic-geometric mean inequality 

c||C/||i2||VC/|U2 <^||C/||2, +«;||vc/||i2, 

we get (|6.5b|) — after dividing through by ||f7||^2. 

Appendix B. Evans function when spectral parameter is large. We follow 
the arguments for such results given in Alexander, Gardner and Jones |I1 p. 186] and 
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Sandstede P^l, p. 24]. We begin by rescaling x — \/\X\x in 
—K, ~K + 1, . . . ,K and taking the limit |A| ^ cx) to get 



for each k 



fei/fe-e'^'s^B-i^Jfe =0. 



Hence we have a decoupled system of 2K + 1 equations of the form 

a. 



02{2K+1) h ® I2K+I 
B ® I2K+I 02(2K+1) 



where B = e^'^^'^^B ^ . This is a constant coefficient linear problem and the so-called 
spatial eigenvalues /i are the roots of the characteristic polynomial 



det{{^J?I ~B)® I2K+1) = {deti^l^I - B)) 



2K+1 



where for any matrix C we have used that det(C Ik) = det(C'"') = (detC) . 
Hence the spatial eigenvalues are zL^/Bn, ±^/B22', each with algebraic multiplicity 
2K + 1. Associated with each of these four basic spatial eigenvalue forms we have the 
{2K + l)-dimensional eigenspaces: 

M = +Vb^- span{e+ (+ : k ^ -K, -K+1,..., +K}, 

M = -Vb^- span{e2j.(-ysI7): k = -K,-K+l, . . .,+K}, 

M = +V^: span{e+^i(+v/S^): k ^ -K, -K + I, . . . , +K}, 

M = -\/b^: span{e2j,+i(-V's^): k = -K,-K -I- 1, . . .,+K}, 

where ei(/3) is the vector with 1 in position i, with (3 in position 2{2K -I- 1) -I- i, and 
zeros elsewhere. Hence as |A| — > cxi, the Evans function has the asymptotic form 

^(^)-d^*U^/^«W B^/-^l2K^, 

= 2det(Bi/2®/2/f+i) 
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